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(<— ^ ' We simulate the Gross-Neveu model in 2+1 dimensions at nonzero baryon 

0^ ■ density (chemical potential /i 7^ 0). It is possible to formulate this model 

with a real action and therefore to perform standard hybrid Monte Carlo 
C^ , simulations with /x 7^ in the functional measure. We compare the phys- 

ical observables from these simulations with simulations using the Glasgow 
Qj I method where the value of // in the functional measure is fixed at a value 

>-^ ■ l^upd- We find that the observables are sensitive to the choice of ^upd- We 

consider the implications of our findings for Glasgow method QCD simula- 

^ , tions at ^ 7^ 0. We demonstrate that the realisation of the Goldstone mecha- 

H I nism in the Gross-Neveu model is fundamentally different from that in QCD. 

We find that this difference explains why there is an unphysical transition 
in QCD simulations at // 7^ associated with the pion mass scale whereas 
the transition in the Gross-Neveu model occurs at a larger mass scale and is 
therefore consistent with theoretical predictions. We note classes of theories 
which are exceptions to the Vafa-Witten theorem which permit the possibility 
of formation of baryon number violating diquark condensates. 



I. INTRODUCTION 

Successful lattice simulations of QCD at nonzero baryon density have yet to be achieved 
and the fundamental obstacle to success is the fact that standard hybrid Monte Carlo tech- 
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niques do not admit the complex functional measure appropriate to the full theory. The 
results of the quenched theory are known to be unphysical [|l]-§]. Attempts have been made 
to study unquenched QCD at // 7^ using the Glasgow method which involves generating 
a statistical ensemble via hybrid Monte Carlo at /x = , reweighting the observables by 
the ratio of the fermion determinant at /i 7^ to that at /x = and performing a Grand 
Canonical Partition Function expansion in the fugacity variable e^^'^ to obtain observables 
for all /i, but the results were not significantly different from the quenched theory. In this 
paper we study the Gross- Neveu model which is one of the few models amenable to hybrid 
Monte Carlo simulations at nonzero density. By implementing the Glasgow method in the 
Gross-Neveu model and comparing observables with those of a standard hybrid Monte Carlo 
simulation we reveal the limitations of the Glasgow method. The results are relevant to QCD 
simulations. Work by Stephanov strongly suggests that any nonzero density simulation 
which incorporates a real path integral measure proportional to det (MM^) is doomed to 
failure due to the formation of a light "baryonic pion" from a quark and a conjugate quark 
P^J^. Since an earher simulation of the Gross-Neveu model at /i 7^ did not exhibit 
such a pathology an explanation is required. We provide a solution to this puzzle later in 
this paper. 

Let us briefly review the current status of QCD simulations at nonzero density. In QCD 
the fermion determinant is complex for chemical potential /x non-zero, therefore generating 
an ensemble by Monte Carlo at /i 7^ is impracticable. The fermion number density, 
Jo = ('?/^7o'?/'), measures the excess of quarks relative to anti-quarks and is expected to start 
to rise from zero at some value, \io. This onset, Ho corresponds to the point where the 
phase of nuclear matter is more energetically favourable than the vacuum state (for which 
Jo = 0). Chiral symmetry restoration occurs at some critical /ic and we expect [^ that 
yUo ~ /ic — "^p/S, where rup is the proton mass. Quenched simulations, on the other hand, 
predict ^o — '^w/S suggesting that in the limit where the bare quark mass m — > chiral 
symmetry is restored for any /x 7^ 0. This result is unphysical because the pion which has 
baryon number zero, should not couple to the baryon chemical potential. 

First attempts to simulate full QCD using the Glasgow method |^ also produced perplex- 
ing results 0. On an 8^ lattice with bare quark mass m^y = 0.01 at /9 = 5.1 we found fio — 0.1 
which differs considerably from the strong coupling analysis [0 prediction /i^ — 0.65 for 
(3 = 5.0. Furthermore the scaling of fio with m was consistent with a Goldstone boson 
controlling the onset. This result has motivated an assessment of the effectiveness of the 
Glasgow method via its implementation in a simpler model. 

Some insight into the effectiveness and the possible problems of the Glasgow method was 
obtained by a study of QCD at inflnite coupling [rT|, where we compared results obtained 



by use of the Glasgow method with the quenched and the nearly exact ones available in 
that limit. Glasgow results for fi ~ 771^/2 were found to reproduce the quenched ones in 
quantitative detail. This rules out plausible physical explanations of the early Glasgow 
onset, which include anomalously low baryonic states or thermal excitations of baryons |l^ 



pi]| . In the present study the Glasgow results will be compared with an analogous of the 
quenched results (Toy model) and with the exact results of [§. 

In section II we describe in detail the lattice formulation of the 3d GN model including a 
discussion of the symmetries and the features of the model which make it amenable to lattice 
simulation by the hybrid Monte Carlo algorithm even at nonzero chemical potential. In 



section III we outline the Glasgow Method for the GN model which makes use of eigenvalue 
symmetry relations to construct a Grand Canonical Partition Function expansion in the 
fugacity variable. The special feature of the Glasgow method is that in principle it permits 
simulations using any fixed value of fi in the functional measure in order to obtain observables 
such as fermion number density for an extensive range of fi values ranging from /x = to 
yU > fic- By comparison with the hybrid Monte Carlo simulations where the fi = fiupd in 
the functional measure matches the /i = Hmeas at which the observable is calculated we 
assess the effectiveness of the Glasgow method within the limits of our statistics. In section 
III we present the results of our Gross-Neveu simulation using the Glasgow Method and 
provide a direct comparison with results of a standard hybrid Monte Carlo simulation. In 
section IV we present results of a Toy Model simulation of the GN model. In the Toy 
Model the /i in the functional measure is different from the n at which the observable is 
calculated. When /z in the functional measure is set to 0, the Toy model is the Gross- 
Neveu equivalent of a quenched approximation. When /i in the functional measure equals 
H at which the observables are calculated, the Toy model is exact. The Toy model results 
therefore give an insight into the relevance of including the chemical potential into the 



dynamics. As discussed in [^ , the Glasgow method should improve over the quenched/Toy 
model because of the reweighting. In the low statistics limit, the Glasgow method should 
reproduce the Toy results, while in the Glasgow method should approach the exact results. 
Hence the interest of considering both Toy and exact results to understand the effectiveness 
of the Glasgow method. In Section V of this paper we discuss the current onset, comparing 
Toy, Glasgow, and HMC results. We will show that these onsets are all controlled by the 
lowest excitation in the pseudoscalar compostite spectrum. In contrast with QCD, this is 
not the Goldstone pion, but a massive state. In section VI of this paper we consider a 
question which is not specific to the lattice simulation method but is highly relevant to 
our understanding of the physics of chiral symmetry restoration. We provide evidence via a 
hybrid Monte Carlo simulation that the realisation of the Goldstone mechanism in the Gross- 
Neveu model (specific to theories with four-fermion interactions) is fundamentally different 
from the familiar mechanism in QCD. We show that the reason for this is that in the Gross- 
Neveu model the dominant contributions to the Goldstone pion come from the disconnected 
diagrams. This fact precludes the existence of a light baryonic pion in /i 7^ simulations 
incorporating four-fermion interactions. Finally in section VII we discuss the relevance of 
our Gross-Neveu model study to QCD simulations at nonzero chemical potential. 

II. LATTICE FORMULATION OF THE 3D GN MODEL 

The fermionic part of the lattice action we have used for the bosonized Gross-Neveu 
model with U(l) chiral symmetry is given by [§ 
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Here, Xi a-nd Xi are complex Grassmann- valued staggered fermion fields defined on the lattice 
sites, the auxiliary scalar and pseudoscalar fields a and vr are defined on the dual lattice 



sites, and the symbol (x, x) denotes the set of 8 dual sites x adjacent to the direct lattice 
site X. N is the number of staggered fermion species and 1/ g'^ is the four-fermi coupling. 
The symbol e{x) denotes the alternating phase (— l)^o+a;i+a:2^ The auxiliary boson fields a 
and TT are weighted in the path integral by an additonal factor corresponding to 
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The fermion kinetic operator }A at non-zero density is given by 
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where m is the bare fermion mass, ^ is the chemical potential, and riy{x) are the Kawamoto- 
Smit phases (_i)^'o+-+2:-'-i_ 

The model can be simulated using the hybrid Monte Carlo algorithm |T^, in which 
complex bosonic pseudofermion fields $ are updated using the action $'''(M^M)~^$. After 
integration over $ the measure in the functional integral is therefore manifestly real and 
given by 
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where we have used the fact the only complex entries of M occur on the diagonal. Therefore 
the simulation includes A^ "white" flavors of fermion with positive axial charge (ie. coupling 
to +171) , and A^ "black" flavors with negative axial charge. If we include the global chiral 
symmetry valid for m ^ 



x(a;) ^-* exp(zQ;£:(x))x(a;) ; xi.^) ^~^ exp(iQ;£:(a;))x(a;) 
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then we see that the model at non-zero lattice spacing has the symmetry: \5 {N)y ®\] {N)v ® 
U(1)a- It is the U(1)a symmetry (|2.5| ) which is broken, either spontaneously by the dynamics 
of the system, or explicitly by a bare fermion mass. In the continuum limit, it is possible 
to recast the model in terms of Nf = 4N flavors of four-component Dirac spinors i/JjI/j 
||14j| , ie. A^//2 white flavors and A^//2 black; the global symmetry group should enlarge 
to U(A^//2)y ® U(A^//2)v ® U(l)^. The continuum action approximated by ( |2.1| ) after 



integration over the auxiliary fields a and vr is | |1 5| | : 
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The relevant features of the 3d Gross-Neveu (GN) model for /i 7^ studies are that it 
has a chiral transition with a massless pion in the broken phase and it can be formulated 
such that the fermion determinant is positive definite for /i 7^ 0. 



III. GLASGOW METHOD FOR THE 3D GN MODEL 

The key feature of the Glasgow Method is that the chemical potential fiupd appearing in 
the functional measure is fixed to a constant value while the associated observables are cal- 
culated at the appropriate fimeas which can be varied at will. The observables are calculated 
from logarithmic derivatives of a Grand Canonical Partition Function (GCPF) expansion 
in the fugacity variable e^^'^ where T is the temperature. To compensate for generating 
the ensemble at fixed fiupd, the observables are "reweighted" by a ratio Rrw of fermion de- 
terminants. This procedure of reweighting also serves to make the method exact i.e. the 
expression for the GCPF from which the observables are obtained is formally correct in the 
limit of infinite statistics. In practice it is necessary to assess the effectiveness of the Glasgow 
Method via direct comparison with a standard hybrid Monte Carlo simulation in which fiupd 
is variable and identical to the fimeas appearing in the expression for the observable. 

Consider the expression for the GCPF in the Gross-Neveu model: 

Z{jj,) = f[da][dTT] det(M(a,7r,/i,m)) e"^""^ (3.1) 

where M is the fermion determinant, Saux contains the auxiliary fields a and n, m is the 
bare/current quark mass and /i is the chemical potential. The GCPF (for fixed m) can be 
rescaled and expressed as an ensemble average of det M at some fixed value /i = ^upd'- 
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where angled brackets denote an average over an ensemble generated with chemical poten- 
tial ^upd, and Hmeas IS the chemical potential for which the physical observables are to be 
calculated. Note that for example generating the ensemble at Hupd = would allow us to 
circumvent the problem of a complex action in a Monte Carlo simulation. 

For optimum efficiency of the Glasgow method we require a large overlap between the 
ensemble generated using det M^fi^pd) in the functional measure and the exact ensemble 
generated using det M{nrneas)- Let us define a reweighting factor Rrw = Jet m'"'"'"rn) ■ '^^^ 
relative magnitude of the factor Rrw configuration by configuration gives a measure of the 
overlap. If there is poor overlap between the simulated ensemble and the true ensemble it 
is conceivable that only a small fraction of the configurations will contribute significantly to 
Z (those where Rrw is large in magnitude) in which case extremely high statistics would be 
required to extract sensible results. 

The full lattice action for the bosonized GN model with U{1) chiral symmetry is given 
in ( p.l|) . The functional measure used in the HMC algorithm is det (MM'''). To formulate 
the Glasgow Method for this model we express the Dirac fermion matrix, M and a related 
matrix M, in terms of matrices G and V where G contains all the spacelike links while V 
(V^) contains the forward (backward) timelike links. Note that for the Gross-Neveu model 
det (MM) = det (MM^). 

2tM,y{fi) = Y,y + G,y + V^yC^ + ^^e"^ ; -2iM.,y{fx) = Y^ + G,y + V^ye^ + r^e-^ (3.3) 

The term describing the Yukawa couplings of scalars to fermions is given (in terms of the 
auxiliary fields a and tt on dual lattice sites x) by 



Y^y = 2t{m + - J2 (ff(5)+«£7r(5)))(5^j/. (3.4) 
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The determinants of these fermion matrices are related to that of the propagator matrix 
P (following Gibbs ^y. 

--(-"-/" o) (-) 

which is a matrix of dimension 2n^nt in this model, for a n"^ x nt lattice. Note that V is 
an overall factor of P. The inverse of the propagator matrix is 

p-'-v'(°_g-y) -d (^-')' = ( _°i _a \ y, ) y (3.6) 

Note that Y^ = Y* 

The determinants of P and M are simply related: 



det(P - e"'^) 



-GV -YV - e-^ V 

-V -le-^' 



= det {GVe-^" + YVe-^" + e'^^ + V^ 

= det {{Ge-^" + Ye"" + V^e-'^" + v) v) 

= e"'^"'"* det (G + r + V^e'^ + Ve^") 

= e"'^"'"* det(2M) (3.7) 

where we have used det V = 1 and V'^V = VV'^ = 1. Similarly the determinants of P and 
M are simply related. 

We shall expand det (MM) as a polynomial in e"'^ and in so doing we make use of two 

symmetries of the eigenvalues. First consider the transformation A = in the space 

of the propagator matrix. One can easily show that 

A(P-)W=(-^^-^' J)=P (3.8) 

Thus if A is an eigenvalue of P then ^ is an eigenvalue of P. The presence of A* and ^ in 
the spectrum then follows from the fact that det (MM) is real for arbitrary fugacity. 
We now perform our GCPF expansion: 

det(MM) = e^"^'""^ det(P - e'^) det(P - e"^) (3.9) 

= e^-^-'^ n {>^^ - e"") [j: - e-^j (a* - e-^) i^- - e"^j (3.10) 

_g2n.2„,^-Q h _ g-M(A. + _) + e-2M h _ e-A'(A* + — ) + e-^M (3.11) 
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The above GCPF expansion incorporates all of the eigenvalue symmetries of the model. 
Provided the (a„) are determined to sufficient accuracy, we can measure the averaged charac- 
teristic polynomial over the ensemble generated at any fixed /i = ^upd using a hybrid Monte 
Carlo algorithm, and use this to provide an analytic continuation §] for the GCPF to any 
non-zero /i. In order to obtain the fugacity expansion we must determine the eigenvalues of 
PP. 

Since the matrix V is an overall factor of P, consider the effect of multiplying a timelike 
link by e^'^*. We can perform a unitary transform to spread this over all timelike links so 
that a new symmetry emerges: 

V — ^ V' X element oiZ{nt) (3.14) 

This is then transferred to the eigenvalues, A 

det(P-Ai) = (3.15) 

Therefore the eigenvalues themselves have a Z{nt) symmetry. This Z{nt) symmetry holds 
configuration by configuration. As a consequence of this symmetry, the characteristic poly- 
nomial for P is a polynomial in e^"' with [An'j. + 1) real coefficients. Thus we obtain an 
expansion for the GCPF in the fugacity, e^/"^ 

Z oc E(«")(e''"* + e"^"')" = Y. e-(^"-"^)/^ (3.16) 

The major computational task in performing the GCPF fugacity expansion is the deter- 
mination of all of the eigenvalues of PP. It is more efficient to diagonalize {PP)^^ than 
to diagonalize PP, therefore we exploit the Z{nt) symmetry of the eigenvalues . This in- 
troduces a Z{nt) degeneracy of the eigenvalues and effectively reduces the dimension of the 
matrix to An^rit to 4n^. Note that although PP is a sparse matrix (PP)"'^ will be dense. The 
expansion coefficients (a^) are evaluated in the simulation and thermodynamic observables 
can be obtained from derivatives of \nZ. In particular the number density is defined by: 



(Jo {fi,ma)) = lim 
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T 91n(Z(/i,m)) 
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and in terms of the GCPF expansion 



(Mil)) = "="'"- . (3.18) 

Since the 3d GN U{1) model has a positive definite functional measure we can choose 
any desired fiupd and investigate the infiuence of this choice on the observables. Standard 



hybrid Monte Carlo (HMC) simulations showed a clear separation of the Goldstone pion 
mass scale and the critical chemical potential {^c) where chiral symmetry was restored i.e. 
Hc» vHt^jI. Thus the lattice Gross- Neveu model is an example of a theory where there is 
no manifestation of a light Goldstone pion carrying non-zero baryon number even when the 
model is simulated with a real path integral measure proportional to det (MM^). 

We set out to investigate whether by applying the Glasgow Method for a given [lu-pd we 
could obtain from the GCPF expansion coefficients a fermion number density as a function 
of [I which would match the results of the standard hybrid Monte Carlo simulations. Of 
particular interest was whether we would observe the discontinuity in the number density 
at yUc ^ w^7r/2 via the Glasgow Method. 

We performed standard hybrid Monte Carlo and Glasgow Method simulations on 16^ 
lattices at a four-fermi coupling of l/(?^ = 0.5 and m = 0.01. To minimise the number 
of GCPF expansion coefficients we simulated the Glasgow Method with A^ = 1 flavours of 
staggered fermion. (Note that the results in [^ are appropriate to iV = 3.) In fig |^ we 
compare number densities from standard hybrid Monte Carlo simulations with A^ = 1 and 
A^ = 3. For N = 3 fic = 0.725(25) while m^/2 = 0.18(1) g]; however for A^ = 1 there 
is a shift such that /ic — 0.6 i.e. fewer fermion flavours implies earlier chiral symmetry 
restoration. 

The simulation results demonstrate that the chemical potential fiupd at which the statis- 
tical ensemble is generated has a strong influence on the thermodynamic observables of the 
simulation. Fig. |^ shows the number density for the standard HMC simulation and for three 
simulations using the Glasgow method: one for fi^pd = 0.0 another for fi^pd = 0.625(~ /ic) 
and flnally with fi^pd = 0.65 (> /ic) The discontinuity at ^c associated with the fermion losing 
dynamical mass, which is apparent in the standard hybrid Monte Carlo data is not consis- 
tently reproduced by the Glasgow algorithm - the results depend strongly on the choice of 
fiupd- In fact the discontinuity associated with the chiral transition can be extracted from 
the GCPF expansion coefficients only when fiupd — fJ^c- 

In this simple model the constraints on the effectiveness of the reweighting must come 
from distribution of the a field. The histograms of measurements of a for the two different 
updates are shown in Fig. |^. For fi^pd = the sigma field is sharply peaked about a mean 
value associated with a large dynamical mass whereas for fiupd = 0.7 the distribution is 
broader although there is still no clear evidence for a two-state signal for a fermion with 
and without dynamical mass. Despite the absence of a two-state signal in the sigma fields 
the reweighting still refiects the discontinuity in the fermion number density although this 
discontinuity is sharper in the standard HMC simulation. 

A. Free Gas behaviour 

We shall compare Jq as a function of /i obtained from the Glasgow method with the 
number density of a lattice gas of free fermions with mass rrif = m+ (a). Firstly we consider 
the simulation with fiupd = 0.0 where we found Jo rises smoothly from ^o to saturation and 
secondly we consider the simulation with fiupd = 0.625 where we saw a discontinuity in Jq 
at a value of fi which was consistent with the fie determined by the standard hybrid Monte 
Carlo where fiupd = A^meas for each value of fi simulated. 



This comparison with free lattice gas number density is useful because it gives us a 
qualitative idea of whether the Jq obtained for a given ^^pd is revealing the transition from 
the chirally broken phase where rrif is large to the chirally symmetric phase with ruf ~ m. In 
fact as we shall see in section V the Toy Model number density is very well approximated by 
a free lattice gas of fermions therefore any significant departure from the free gas behaviour 
for the Glasgow Method number density must be due to the reweighting factor Rrw 

The fermion number density for a free Fermi gas is obtained by summing over the finite 
set of Matsubara frequencies |jl8| that arises in the lattice thermodynamics. The appropriate 
expression for the number density, Jq in terms of the momentum sums for a 3-dimensional 
free lattice gas is :- 



•JQfree — ^^ / 7^34 :^=;2 — „;„2 „ — , „;„/„ — rTTTVTTU 1,0. iyj 



(Pp sin(pi + ifi) cos{pt + ifi) 
-TV 2iT^ T4=i sin^ Pi + sm{pt + ^/i) + mj 

where pt is associated with the temporal direction. 

Recall that (a) is a measure of the physical fermion mass. The results of our Jq compar- 
ison are plotted in Fig. ^. In our GN simulations using the Glasgow method for fiupd = 0.0 
the fermion number density was consistent with a gas of free fermions with a dynamical mass 
rrif = m+ (cr)|^=o- As we shall describe in section V, this is the same behaviour we observed 
for the Toy model. However for fiupd = 0.625 there is a clear indication of a discontinuity 
at fic and for n > ^c the results are consistent with a free gas of fermions with dynamical 
mass mfe + (o")|^=o.625- We found that for n^pd = 0.65 i.e. as the chemical potential in the 
functional measure increased beyond /Zc the discontinuity disappeared again. This suggests 
that within the limits of the statistics of a typical lattice simulation, the inclusion of the 
reweighting term Rrw in the expression for the GCPF is insufficient to project onto the 
ensemble appropriate to iimeas except when /Xmeas — ^^upd■ 

B. Lee- Yang Zeros 

We have seen that the Glasgow Method fermion number density only gives a reliable 
estimate for /Xc from the number density when ^^pd — f^c'- however we also expect a signal 
for criticality from the partition function zeros. According to the theorems of Lee and 
Yang ||19|, the phase transitions of a system are controlled by the distribution of roots 
of the GCPF. A phase transition occurs whenever a root approaches the real axis in the 
infinite volume limit. In practical lattice simulations we are certainly not close to the 
thermodynamic limit, however as the lattice volume is increased we expect the zero with 
the smallest imaginary part to approach the real axis, converging to the critical value of 
the appropriate parameter. The Lee- Yang zeros in the complex ^ plane are the zeros of 
eqn. ( |3.1| ) and their distribution should refiect fic- The zeros in the complex /i plane for 
the fiupd = 0.0, 0.5, 0.55, 0.625, 0.65, 0.85 are plotted in Fig. |[ The first point to note is that 
the zeros form a distinctive and elaborate pattern. Even for this simple model it would be 
difficult to predict analytically the distribution of zeros however some very clear features 
emerge when we compare the distributions appropriate to different fiupd- For fi^pd < yUc i-e. 
f^upd = 0.0, 0.5, 0.55 there is an arc of zeros which intersects the real axis at /i ~ 0.6 ~ fic- 
For these 3 values of fiupd we saw no evidence of a discontinuity in Jq at // = 0.6. Now 



at fiupd = 0.625 there is clear change in the pattern of zeros. What was an arc of zeros at 
/i ~ 0.6 now becomes a broader vertical line and two additional lines of zeros sit outside 
the main body of the distribution at /i ~ 0.5. Recall that for fiupd = 0.625 we did see 
a discontinuity at /ic in Jo- As we increase fiupd beyond fie i.e. for fiupd = 0.65, 0.85 the 
vertical line of zeros migrates to higher fi and the main zeros distribution extends further 
into the low fi region of the complex fi plane. For fiupd = 0.85 the vertical line of zeros sits 
at yU ^ 0.85 i.e. far from fie- 

IV. TOY MODEL VERSUS STANDARD HYBRID MONTE CARLO RESULTS 

In this section we discuss an approximate model, which we will refer to as the 'Toy' 
model, which will turn out to give useful insight into the results of Glasgow reweighting 
applied to the GN model. In the Toy Model the chemical potential used in the measure 
/^meos can be difefrent from that used in the update fiupd- It can be seen as the equivalent of 
the quenched approximation when fiupd = 0, and is the exact model when fimeas = fJ'upd- It 
is essentially the Glasgow Method with Rrw = 1. When fiupd = and the statistics is low, 
the Glasgow method results coincide with thos of the Toy model, and should approach the 
exact results in the large statistics limit |TI|]. 

Let us begin by reviewing what is known about the GN model in the mean field approx- 
imation. 

In the mean field approximation, which turns out to be the leading order of an ex- 
pansion of powers of 1/iV/, we solve for the expectation value of the auxiliary scalar (a) 
self-consistently using the gap equation, which reads (with bare mass m set to zero): 



(^) = -g'im = ytrS^m) = l^rg'-^^J^- (4.1) 

Eqn. (|4.1|) can be solved using a lattice regularisation [^ and has a non-trivial solution 
(a) 7^ for gf^ > ^f^ ~ 1.0. The dynamical fermion mass rrif = (a) then defines the model's 
physical scale, and there is a continuum limit corresponding to a continuous chiral symmetry 
restoring phase transition as {a) \0, g^ \ g^. 

The gap equation can also be generalised to /i 7^ 0: 

W(/i) = -9'{^^){f^) = ytrSpifi, (a)(/i)); (4.2) 



once again the solution can be found, either in the continuum [^, or on a lattice |2^ 0. 



The salient feature of the solution is that (o")(/i) remains unchanged from its zero-density 
value (cr)(0) as fi is increased, up to a critical value 

/x,= (a)(0), (4.3) 

whereupon it falls immediately to zero, signifying that chiral symmetry is restored via a first 
order transition. Introduction of a small bare fermion mass m causes a slight softening of 



the transition. The hybrid Monte Carlo simulation results of refs. ||22| [g] confirm that this 
picture is qualitatively correct. 
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Now consider a 3d GN model in which a value of the chemical potential ^upd is used to 
generate the ensemble, and a value fimeas to measure expectation values, with fiupd 7^ f^meas in 
general. We will refer to this as the ' Toy modeV Notice that fiupd = fJ-meas defines the standard 
hybrid Monte Carlo calculation whereas fi^pd = is the equivalent of a quenched study (we 
could indeed refer to the case fiupd 7^ fJ'meas as a 'generalised quenched approximation'). 

In a mean field treatment of the Toy model the a field is constrained to its fi = fiupd 
value, so therefore we simply substitute (a) = {(j) (fiupd) in the fermion propagator, and 
calculate (ipip)- Therefore the Toy counterpart of eqn. ([4.2| ) is 



{o-){fiupd) 7^ -5'^ <i)ll)> iflmeas) = T^tr^F (yUmeas, {ct) (flupd)) (4.4) 

In the mean field approximation, the above considerations suggest that all Toy models are 
equivalent to free fermion models with ruf = {cr{fiupd))- This is confirmed by the numerical 
simulations described below. It is clear that Glasgow method reweighting cannot work in 
the GN model because the ensemble generated with fj,upd < f^c has zero overlap with that 
generated using fi^pd > l^c- This means that all Toy models with n^pd < f^c should be 
equivalent, and should approximate the standard HMC results for fimeas < f^c- However, for 
fimeas > yUc the propagator 5*^ remains that of a free fermion with mj = {(j){fiupd) i-e. the 
fermion appears to retain its dynamical mass therefore we cannot expect a chiral symmetry 
restoring transition in the Toy model. Consider next the Toy model with Hupd > i^c- The 
propagator Sp now always describes a massless fermion, therefore in this case we expect the 
Toy model to reproduce the exact results only for /Xmeas > /^c- 

It is also possible, and instructive, to test these ideas once quantum fluctuations of a 
are included. Accordingly we have performed simulations of the Toy model for ^upd = 0, 
fimeas arbitrary, using the hybrid Monte Carlo algorithm. Fig. ^ compares the results of 
these simulations with the exact simulation results 0, obtained with the same parameters 
(coupling l/g"^ = 0.5, bare mass m = 0.01, lattice size L = 16, A^ = 3 corresponding to Nf = 



12 continuum flavours), and the mean field predictions of (f4.4|) . Both the condensate (V''0) 
and number density (Jo) are plotted as functions oi fimeas- We see that the Toy model with 
fiupd = is well described by a massive free field with m/ = (o")(0) for fimeas < f^c — 0.725, 
and that exact and Toy results agree well in this range; for fimeas > fie, however, while the 
exact results for (tptp) fall steeply, the Toy results exhibit a smooth crossover, approaching 
zero only for fimeas — 1-4. For fimeas > fie the full model is better described by a massless 
free field, corresponding to fiupd > fie, {o') (fiupd) = 0, indicating the restoration of chiral 
symmetry. In all cases the mean field prediction (|4.4| ) reproduces the Toy simulation results 
well, indicating that quantum fluctuations are small due to the relatively large value of Nf. 

We deduce that fie is not observable in the Toy model, since there is no chiral symmetry 
restoration. Note however that the value /Iq at which the onset of (Jq) occurs is the same 
for both Toy and exact models, reflecting the fact that the onset on a finite lattice (and 
perhaps even in the thermodynamic limit) is not associated with the chiral transition. 

We emphasize that this simple analysis works well because at l/g"^ = 0.5 and Nf = 12 
the fluctuations in a are highly suppressed. We verified that the microscopic distribution of 
the a field is very sharply peaked about its constant value. This condition serves to makes 
the model well approximated by the appropriate free fermion field in each phase. Because 
the two phases for fi < fic and fi > fie are clearly different, we anticipate that Glasgow 
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reweighting using Hupd = will not be able to capture the chiral phase transition because 
the overlap between a distributuions from each phase is likely to be very small. The overlap 
would be exactly zero in the infinite Nf limit where fiuctuations are completely suppressed. 
It would be interesting to repeat this exercise closer to the critical point g^ — ^ although at 
this coupling the chiral transition is more difficult to identify. 

V. PSEUDOSCALAR MASS AND REALISATION OF THE GOLDSTONE 

MECHANISM 

Although we have shown that the Glasgow method and the Toy model simulations differ 
in their quantitative predictions for the number density Jo{fi) from the standard hybrid 
Monte Carlo simulations we observe that the onset, Ho in Jq (ie. the value of /i where Jq 
starts to increase from zero which is usually interpreted as the point of chiral symmetry 
restoration) is correctly predicted by both Toy model and the Glasgow method. This is true 
even for simulations with modest statistics. 

An important result reported in [^ was the numerical proof that fio ^ rrivr/S in standard 
hybrid Monte Carlo simulations of the GN model. Thus we had an example of a theory which 
could be simulated at /i 7^ with dynamical fermions and where there were no pathologies 
associated with the fact that there was a Goldstone pion in the spectrum. In the Gross- 
Neveu theory even in the Toy model simulations (which are partially quenched) fio S> m^/2 
in agreement with the full theory. Why then are there pathologies associated with m^ in 
QCD? We set out to gain insight into this fundamental difference between the two lattice 
models. 

We shall discover that the onset in Jq always occurs at fio = mps/^, where mps is the 
mass of the pseudoscalar state measured via GG^ where G^G'^) is the fermion (anti-fermion) 
propagator . In QCD mps is the Goldstone pion, however in four fermion models the 
Goldstone mechanism is realized by the auxiliary field tt and mps is a much heavier state 
which appears to be close to twice the dynamical fermion mass. 

Let us first consider the realisation of the Goldstone mechanism in the Gross-Neveu 
model. In particular we measure and compare the connected and disconnected contributions 
to the Goldstone pion. If the pion has a dominant connected contribution we expect a pole 
in the pseudoscalar propagator corresponding to a light particle with massoc y/m. If, on the 
other hand, the pion has a dominant disconnected contribution the mass scale corresponding 
to the pole in the pseudoscalar propagator will be ~ 2m/ (yU = 0) >> m^. 

Consider the lattice Ward identity for the chiral condensate: 

J2('^^5^iy)^^5^i^)) = H{'tr{Gl^{x,y)G+^{x,y))) - {{tr-i^G{x,x)){tr-i^G{y,y))) 
V y 

(V^(x)V^(x)) 



rriq 



(5.1) 



The pion susceptibility Eqn. (|5.1 | ) has contributions from a connected channel (1st term) 
and a disconnected channel (2nd term). 

Consider the Dirac kinetic operator Ip for staggered fermions. At /x = the relation 
Ip"^ = —p holds. For /i 7^ the fact that e^ (e~^) multiplies the forward (backward) 
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timelike gauge links in the lattice action means that p'^ ^ ~Pi because in the matrix p'^ , 
g-M (^gM^ now multiplies the forward (backward) links. Consequently the propagators G 
(with G = M~^) and G^ are no longer trivially related as they are at /i = 0. Thus the 
pseudoscalar propagator, Gps, is defined by 

Gpsit)=J2G+,ix,t)Gl^ix,t) (5.2) 

X 

Let Xps be the smallest eigenvalue of the Gibbs propagator matrix P (eqn. |3.5|) . It can be 
shown that Xps is associated with the mass pole in Gps and that when we consider a single 
lattice configuration Xps directly corresponds to a Lee- Yang zero and therefore induces a 
singularity in Jo i-e. triggers the rise from zero of the fermion number density. To see this 
observe that det M can be simply expressed in terms of the eigenvalues of P and Z can be 
similarly expressed in terms of its zeros a^ in the e^ plane: 

det M oc n(e'' - Ai) ; Z oc Yl{e^ - ai) . (5.3) 

i i 

Since Z = (det(M)) we see that on a single configuration a^ = Aj. Note however that 
the ensemble averaged a^'s are not in general the same as the ensemble averages of the Aj's. 

The number density on a single configuration Jq ~ 9 In det M/d^ whereas the ensemble 
average is given by Jq ~ 51n(det M)/c?/i. Clearly the chemical potential where the number 
density begins to rise from zero on a single configuration will be determined by the numerical 
value of Xps (i.e. fio = Xps). 



It has been proved by Gibbs [T^ that if we assume that the pion has a dominant connected 
contribution it follows that Xps ~ e~^. In this case we can only envisage achieving a 
physically meaningful result {fio > ^) if «« 7^ Aj. However in the exceptional circumstance 
where the disconnected contribution to the pion is dominant we expect that Ap^ ~ e"™-'' so 
that configuration by configuration fio = fnf > ^. In this section we shall prove that in 
the 3d GN U{1) model the Goldstone pole forms in the disconnected channel therefore the 
state described by Gps no longer corresponds to the Goldstone pion. Instead we find: 

Gps ~ e-2'^/(^=°)* (5.4) 

Since \Xps\ will now correspond to the dynamical fermion mass rrif rather than 171^^/2 we 
have no reason to expect an early onset in the Gross-Neveu model despite the existence of 
a light Goldstone pion in the spectrum. 

In Fig. 1^ we have plotted the propagator for the pion (measured from the auxiliary field), 
the square of the fermion propagator and the pseudoscalar propagator. The results show that 
the pseudoscalar state is very close to twice the physical fermion mass and mps > 2m f >> 
TTij^. What is remarkable is the fact that the auxiliary field pion is considerably lighter than 
mps- This can be understood in terms of the mechanism of chiral symmetry restoration 
in 4-fermion interaction models whereby the disconnected contributions are responsible for 
making the pion light. We will provide numerical evidence for this below. If we measure 
the pseudoscalar propagator in the standard way from GG^ we include only the connected 
contributions to the particle mass. 

Fig. ^ shows the disconnected contribution to the Ward identity for a standard hybrid 
Monte Carlo simulation with A^ = 3 at zero chemical potential. For this simulation we 
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found that the signal for chiral condensate from the noisy estimator was relatively stable as 
shown in Fig. ^ giving {iljijj)/mq ~ 40, therefore we require that the sum of the connected 
and disconnected diagrams be of similar order so that Eqn. ^?1] is satisfied. The connected 
contribution was stable at a value of around 0.5 therefore the dominant contribution must 
come from the disconnected diagram. We found that the disconnected contribution was 
very noisy with large downward peaks. The data suggests that considering the connected 
contribution alone will never be sufficient to satisfy Eqn. p.l| . 

We repeated our measurements for a non-zero chemical potential. We chose /i = 0.5 
thus ensuring that we were still in the phase of broken chiral symmetry. The results were 
consistent with those at zero chemical potential as one would expect. 

As a consistency check observe that we expect the following relationships to hold from 
the equations of motion: 

^^ = _a ; ^V^e^(a;)V55^(y) = — — ^7r(x) (5.5) 

r y 9 y \x J 

The first relationship has been checked and is satisfied by the simulation data. The second 
relation can also be verified. As seen in Fig ^ the characteristic peaks in the disconnected 
contribution are clearly correlated with the square of the pion field. We note that in the 
exact HMC study 0] the pion mass was obtained from the correlator of the auxiliary field 

71. 

The disconnected contribution to the pion susceptibility gives a very noisy signal which 
suggests that a very long run would be required to equilibrate sufficiently to satisfy the 
lattice Ward identity. 

VI. RELEVANCE OF GN STUDY TO QCD SIMULATIONS 

The Glasgow Method for QCD is very similar to the method described in section III 
of this paper however because the QCD action is complex for (/i ^ 0)we are restricted to 
generating the ensemble at fiupd = 0. Furthermore the relationship det MM = det MM^ 
which we used in the lattice Gross-Neveu model does not hold for QCD at /i 7^ 0. We have 
seen in section IV that the GN model observables are sensitive to the choice of fiupd and 
the most accurate results were obtained for fiupd — f^c- Clearly it is impossible to simulate 
QCD at ^upd — fnp/3 therefore we should be aware that very high statistics are likely to 
be required to obtain accurate observables from the Glasgow Method for QCD. It is still 
possible that the Lee- Yang zeros distribution will give a hint for /Xc even with moderate 
statistics and fiupd = 0. Since the Gross-Neveu model is a purely fermionic theory it is 
impossible to assess from our study what infiuence the gauge fields have in QCD dynamics 
in the context of Glasgow method reweighting [|r^ . 

In section V of this paper we have proved that the disconnected contribution to the 
Goldstone pion is dominant in the Gross-Neveu model. In QCD on the other hand we know 
that the pion has a dominant connected contribution which implies that 



Gps{t) oc e 



-rriTrt 



c.f. eqn. p.4|) and this explains the origins of the early onset of the chiral transition in 



QCD associated with a light baryonic pion. The existence of the baryonic pion in the 
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Gross-Neveu model does not lead to an early onset of the cliiral transition because it is 
the disconnected contributions which make the pion light i.e. the Goldstone mechanism is 
realised in a fundamentally different way. 

VII. CONCLUDING REMARKS 

In this paper we have found an example of a model where the inclusion of the chemical 
potential into the dynamics is essential to obtain exact results. 

We have proved that the physical observables of the Glasgow Method are sensitive to the 
choice oi fiupd- The optimal choice for revealing critical behaviour is fiupd — fJ^c, a choice which 
is impractical for QCD simulations. In QCD the pseudoscalar channel pole is formed from 
connected diagrams corresponding to G^^{t)G_^{t). A fermion pair with nonzero baryon 
charge (baryonic pion [^]) forms from a quark and a conjugate quark (i.e. a fermion 
associated with M^ but with the sign of /x reversed) and as a consequence the mass scale of 
the lightest baryon {mn,) is mT,/2 not ra^j?). 

This induces the unphysical early onset of chiral symmetry restoration in quenched QCD. 
In the full QCD it is conceivable that the phase of the determinant will eliminate the equality 
between quarks and conjugate quarks to allow us to recover the physical result m/^ = mp/3. 
In principle the factor R^^ in the Glasgow algorithm should influence the observables to 
take account of the phase of the determinant. The Glasgow method QCD results [§] still 
suggest m/fe ~ m-,j/2. This is likely to be explained by ineffectiveness of the reweighting. 
The reweighting in the Gross-Neveu model has been shown to be ineffective and this is 
not surprising because in this case the statistical ensembles characterising the two phases 
are non-overlapping. The Glasgow method may however be more effective in models more 
sophisticated than Gross-Neveu and in other dynamical regimes such as at high temperature. 

In 3d GN U{1) the Goldstone mechanism is realised by a pseudoscalar channel pole 
formed from disconnected diagrams and the state G+^{t)G_^{t) yields a bound state of mass 
2nif{fi = 0) which is considerably heavier than the pion. As a consequence even when we 
consider individual configurations in this model we find fio — fJ'c^ ni.,^/2. 

The different realisations of chiral symmetry breaking in the GN model and QCD have a 
simple physical origin. For QCD, a vector-like theory, like charges repel, so that light states 
only form between qq pairs (the pathologies in simulations of QCD with /x 7^ occur because 
of the influence of conjugate quarks in the measure; in this case light states may be formed 
from qq"^ pairs). In theories with Yukawa- like interactions such as the GN model, however, 
in the Born approximation all interactions are attractive, and thus one might expect light 
states made of qq, qq and qq. Only the qq system has contributions from the disconnected 
diagram, however. As we have seen, binding in the connected channel is insufficient to make 
light states. Therefore we see a natural relation between the dominance of the disconnected 
diagram and the absence of light states made of two quarks in the spectrum. 

It is also worth reconsidering why simulations of the GN model seem not to be afflicted 
with the problems observed in QCD. As we saw in Section ^, the price of simulating with 
a real measure det (M'^'M) is that the model contains both "white" and "black" fermions 
with opposite axial charges. We might therefore worry about attractive interactions between 
white and black fermions and the possibility of light bound state in the G^'"^'^G''^°f^ channel 
The reason that no spurious onset at /i = m^/2 occurs, and that the GN simulations yield 
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results in good agreement with the mean field predictions, is intimately related to the fact 
that no light state occurs in this channel indeed the white-black interaction due to pion 
exchange in this model has the opposite sign and is thus repulsive: this is in contrast to 
the situation in gauge theories with real measure, where the interaction between quarks and 
conjugate quarks is attractive. 

There are interesting parallels with the standard discussion of the Vafa-Witten theorem 
PB| , forbidding the spontaneous breaking of global vectorlike symmetries such as isospin or 
baryon number in field theories with a positive definite measure in the path integral (this 
includes gauge theories with zero chemical potential but not Yukawa or GN models, for which 
the measure is in general complex). Vafa and Witten do in fact discuss a 'white/black' model 
with a real measure similar to ours, but with Yukawa couplings to a pseudoscalar field only. 
Their analysis does apply in this case, and a fiavor- violating (^ip^hite^biack^ condensate is 
forbidden, the reason being simply that the interaction between white and black particles is 
repulsive. 

Once Yukawa couplings to scalar degrees of freedom are introduced, however, as in the 
GN model, the inequalities necessary to prove the Vafa-Witten theorem no longer hold. 
There appears to be no fundamental obstruction, therefore, to the generation of a baryon- 
number violating diquark condensate l^^p'^hite^white^ -^^ ^^^^ model, despite its measure being 
real; it has been recently suggested that such a condensate forms in QCD at high density 



24 1. A lattice study of this phenomenon is in progress [25 



VIII. ACKNOWLEDGEMENTS 

We thank J.F. Lagae for his comments related to measurement of the pion mass. MPL, 
SEM and SJH received support from EU TMR contract no. ERBFMRXCT97-0122. JBK 
thanks the National Science Foundation ( NSF-PHY96-05199 ) for partial support. 



16 



REFERENCES 

[1] I. Barbour et.a/.,Nucl. Phys. B275, 296 (1986). 

[2] C.T.H. Davies and E.G. Klepfish Phys. Lett. 256B, 68 (1991). 

[3] M. Stephanov, Phys. Rev. Lett. 76, 4472 (1996). 

[4] A. Gocksch, Phys. Rev. D37 (1988) 1014. 

[5] M.-P. Lombardo, J.B. Kogut, D.K. Sinclair, Phys. Rev. D54,2303 (1996). 

[6] S. Hands, S. Kim, J.B. Kogut, Nucl. Phys. B442, 364 (1995). 

[7] M.A. Halasz, A.D. Jackson, R.E. Shrock, M.A. Stephanov, J.J.M. Verbaarschot, [hep- 



[9 
[10, 

[11 

[12 
[13 

[14 
[15: 
[16 
[17 
[18 
[19 
[20 
[21 
[22 
[23 
[24" 



ph/9804290 



LM. Barbour, A.J. Bell, Nucl. Phys. B372, 385 (1992). 

I.M. Barbour, J.B. Kogut, S.E. Morrison, Nucl. Phys. B (Proc. Suppl.) B53, 456 (1997). 

N. Bihc, K. Demeterfi, B. Petersson, Nucl. Phys. B377 615 (1992). 

I. Barbour, S.E. Morrison, E. Klepfish, J.B. Kogut, M.P. Lombardo, Phys. Rev. 

D56,7063 (1997). 

A. Hasenfratz and D. Toussaint, Nucl. Phys. B371 (1992) 539 

S. Duane, A.D. Kennedy, B.J. Pendleton and D. Roweth, Phys. Lett. B195, 216 (1987). 

C.J. Burden and A.N. Burkitt, Europhys. Lett. 3, 545 (1987) 

S.J. Hands and J.B. Kogut, |hep-lat/97050T5i (1997). 



P. Gibbs, Phys. Lett. B172, 53 (1986); P. Gibbs, Phys. Lett. B182, 369 (1986). 

Susan E. Morrison, Nucl. Phys. A642, 269 (1998). 

H. Matsuoka, M. Stone, Phys. Lett. 136B, 204 (1984). 

T.D. Lee, C.N. Yang, Phys. Rev. 87, 404 (1952); Phys. Rev. 87, 410 (1952). 

S. Hands, A. Kocic, J.B. Kogut, Ann. Phys. 224, (1993) 29. 

B. Rosenstein, B.J. Warr and S.H. Park, Phys. Rev. D39, 3088 (1989). 
S.J. Hands, A. Kocic and J.B. Kogut, Nucl. Phys. B390, 355 (1993). 

C. Vafa and E. Witten, Nucl. Phys. B234, 173 (1984). 

M. Alford, K. Rajagopal, F. Wilczek, Phys. Lett. B422, 247 (1998); R. Rapp, T. 
Schaefer, E.V. Shuryak, M. Velkovsky, Phys. Rev. Lett. 81, 53, (1998). 

[25] S.J. Hands and S.E. Morrison, |hep-lat/9807"033 . 



17 



FIGURES 



0.50 



• TOY simulation 
X exact HMC simulation 
lattice free gas m=m^+<a>l i,,, 




0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 1. 
CHEMICAL POTENTIAL (\i) 



2.0 



I.U 




^-x ...^F""^ 






• TOY simulation /•* / 
X exact HMC simulation / / 








- 


0.8 


free lattice gas m=m^ / L 


o 

3 




/ / 

h 1 

1 i 

1 I 

h I 


DENSITY 

o 


- 


DC 

LU 0.4 

CD 

■z. 


" 


.^ 1 
/- 1 


0.2 
nn 




/ / 

^ X V 

BTl— .n.r^ L 4B » 



0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 
CHEMICAL POTENTIAL (^) 



FIG. 1. < ipip > (upper) and Jq in the full and Toy GN model at /3 = 0.5, niq = 0.01, L = 16, 
A^ = 3. The data for the full model fi < 0.9 are from [6] 
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FIG. 3. Fermion number densities comparing simulations for three different values of Hupd with 
the exact hybrid Monte-Carlo simulation on a 16^ lattice with l/g"^ = 0.5, iV = 1 and m = 0.01. 
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FIG. 8. Measurements at ^ = 0.0 showing the disconnected contributions to the pion sus- 
ceptibihty (upper) on a 16^ lattice with 1/g'^ = 0.5, N = 3 and m = 0.01 which we expect to be 
correlated with the auxiliary vr field via an equation of motion. 
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FIG. 9. Connected contribution to the pion susceptibility (upper) and the chiral condensate 
(lower) at ;u = 0.0 on 16^ lattice with l/g"^ = 0.5, iV = 3 and m = 0.01. 
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FIG. 10. Measurements at ^ = 0.5 showing the relative magnitudes of the disconnected (top) 
and the connected (bottom) contribution to the pion on a 16^ lattice with l/g^ = 0.5, A^ = 3 and 
m = 0.01. 
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